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Abstract. A Green’s function method for obtaining an estimate of the ocean circu- 
lation using both a general circulation model and altimetric data is demonstrated. 
The fundamental assumption is that the model is so accurate that the differences 
between the observations and the model-estimated fields obey a linear dynamics. 
In the present case, the calculations are demonstrated for model/data differences 
occurring on very a large scale, where the linearization hypothesis appears to be 
a good one. A semi-automatic linearization of the Bryan/Cox general circulation 
model is effected by calculating the model response to a series of isolated (in both 
space and time) geostrophically balanced vortices. These resulting impulse re- 
sponses or “Green’s functions” then provide the kernels for a linear inverse problem. 
The method is first demonstrated with a set of “twin experiments” and then with 
real data spanning the entire model domain and a year of TOPEX/POSEIDON 
observations. Our present focus is on the estimate of the time-mean and annual 
cycle of the model. Residuals of the inversion/assimilation are largest in the western 
tropical Pacific, and are believed to reflect primarily geoid error. Vertical resolution 
diminishes with depth with 1 year of data. The model mean is modified such that 
the subtropical gyre is weakened by about 1 cm/s and the center of the gyre shifted 
southward by about 10°. Corrections to the flow field at the annual cycle suggest 
that the dynamical response is weak except in the tropics, where the estimated 
seasonal cycle of the low-latitude current system is of the order of 2 cm/s. The 
underestimation of observed fluctuations can be related to the inversion on the 
coarse spatial grid, which does not permit full resolution of the tropical physics. The 
methodology is easily extended to higher resolution, to use of spatially correlated 
errors, and to other data types. 


1. Introduction 

The TOPEX/POSEIDON (T/P) altimetric satellite 
is providing oceanographers for the first time with a 
continuing global database for describing and under- 
standing the large-scale ocean circulation. The satellite 
was launched in August 1992 and ever since has been 
measuring the global sea surface height every 10 days 
with an unprecedented accuracy using a radar altimeter 
system along repeating ground tracks (see the Journal 
of Geophysical Research special issues on T/P, volume 
99, number C12 1994, and volume 100, number C12, 
1995). With this data set, and with other global-scale in 
situ data from programs such as the World Ocean Cir- 
culation Experiment, one finally can begin describing 
the ocean circulation quantitatively on a nearly day-to- 
day basis, rather than as a vague climatological mean. 
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There is another major reservoir of knowledge of the 
ocean circulation: the general circulation models, which 
attempt to summarize understanding of ocean physics 
in the form of numerical codes. The problem we address 
here is to find a practical means for combining the in- 
formation from diverse observations with the dynamics 
inherent in ocean general circulation models (OGCMs) 
so as to produce best estimates of the oceanic state at 
any given time and place accounting properly for errors 
in both models and observations. This problem is a 
challenging one because model and observational errors 
contain complex space/time structures and the number 
of variables required to describe the oceanic state at any 
one time (the “state vector”) is enormous. The general 
problem we are addressing is one of estimation theory 
(meteorologists would call it “assimilation”). 

Algebraically, an OGCM can be written in canonical 
form as 

x(t + At) = r(x(t),q(t),t), (1) 

where x(t) is the state vector at discrete time f, T rep- 
resents the operator stepping the model forward in time 
starting from a prescribed initial condition x(to), and 
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q(t) represents externally specified boundary conditions 
and sources and sinks. 

Most oceanographic measurements, including alti- 
metric ones, are at least approximately a linear com- 
bination of the model state vector, e.g., velocity, tem- 
perature, and salinity, but are contaminated by noise. 
Such measurements can be written as 

y(t) = E(t)x(t) + n(t), (2) 

where E is the “observation matrix” relating the model 
state vector to observables. Normally, E is very sparse 
because observables usually involve only local subsets or 
related elements of x(t). Altimetric surface height ob- 
servations C, in particular, are obtained only at the sea 
surface, measuring the dynamically induced surface el- 
evation (or the equivalent surface pressure p s = gp() 
along subsatellite tracks sequentially in time, while 
floats or tomographic experiments provide subsurface 
data along float trajectories or tomographic ray paths. 
The precise structure of E and the existence of its in- 
verse is a very important question which, in the context 
of altimetry, is directly related to the general issue of 
the observability of the ocean circulation through sea 
surface height measurements [ Fukumori et af., 1993]. 
In the terminology of inverse problems, observability is 
equivalent to the state estimation problem having full 
rank. 

The problem we wish to solve is to find an estimate 
x(t) of the state vector of the ocean subject to model 
dynamics, and its uncertainty P(t), given observations 
y(t), with noise covariance R(t), and a model (1) with 
uncertainty Q (<). Most estimation methods addressing 
this problem require the minimization of a quadratic 
function measuring the model minus data misfit, 

J = £>(«) - Ex(*)) T R(*) -1 (y(t) - Ex(i)) (3) 

t 

subject to the model physics constraints. (Here x(i) 
represents the estimate of the state vector that results 
from combining the model and observations.) There are 
many schemes for minimizing J subject to the model in 
(1), through either constrained or unconstrained opti- 
mization methods. Optimization schemes known to be 
useful, however, overwhelm present-day computer facil- 
ities for realistic problems on basin scales and larger. 
For this reason, many suboptimal approaches are being 
studied, including state reduction, the use of steadys- 
tate sequential algorithms, and the like [e.g., Fukumori 
et a/., 1993; Marotzke and Wunsch , 1993]. 

In addition to the problem of state vector dimension, 
the oceanographic case is complicated by the nonlineari- 
ties in (1). To deal with both the size of the problem and 
the nonlinearity, we use here a method which is based 
upon a form of Green’s functions for a linearized model 
associated with the OGCM. The fundamental assump- 
tion is that modern nonlinear OGCMs driven by ob- 
served meteorological fields will be sufficiently accurate 
in simulating the ocean that the discrepancies between 
the computed and observed oceanic state obey linear 


physics. Almost all optimization methods which are 
practical with large systems (filters, smoothers, Pon- 
tryagin principle/adjoints) are based upon linearization 
of the model dynamics about a reference state. Finding 
the linearization can be an onerous job; one of our goals 
here is to render this job less demanding of programmer 
skill than it normally is. 

As discussed in more detail below, the OGCM is used 
twice: once to compute, using the best available initial 
and boundary conditions, a model first estimate of the 
oceanic state. In a separate calculation, the model first 
estimate is then empirically linearized by perturbing it 
with a succession of simple disturbances on a coarser 
grid. This new, coarse resolution, linear model is used 
to represent the differences between the large-scale data 
and the original full OGCM during the estimation pro- 
cedure. The result is a linear inverse problem for a set of 
disturbance coefficients, which is then solved. The orig- 
inal full nonlinear model can be iteratively improved 
if necessary. This linearization procedure, if it can 
be shown to be effective, would serve several practical 
purposes, including separating the problem of making 
model improvements in the OGCM from those of esti- 
mation, and using the numerical power of a computer 
to produce the model linearization, which is otherwise 
a major chore. 

Note that this method is related to one used previ- 
ously by Wunsch [1988] and Memery and Wunsch [1990] 
for estimating boundary conditions using transient trac- 
ers. Those authors determined the tracer boundary 
conditions which best fit observed tracer distributions 
(in space and time) in the interior ocean. However, as 
will become clear later, the tracer approach is comple- 
mentary to ours: we aim to estimate the interior ocean 
state which best fits observed, but noisy, boundary con- 
ditions. 

There is no guarantee that the linearization is valid: 
the OGCM/data differences must be assessed in detail. 
Furthermore, to usefully combine observations with a 
model, one must be able to show that their differences 
are consistent within the uncertainty estimates for both. 
A full model/data consistency check is an important 
study in its own right and will be presented elsewhere 
[e.g., Stammer et a/., 1996]. 

Here we seek only to make it plausible that the de- 
viations of data and model satisfy a linear set of equa- 
tions. Stammer and Wunsch [1994] show the differences 
between a 3-year mean sea surface height field obtained 
from the global Semtner and Chervin [1992] model with 
1/4° horizontal resolution relative to a 1-year mean of 
T/P sea surface height observations. Deviations are 
generally found in the 10 cm range on large spatial 
scales with length and velocity scales of the order of 
500 km and 10 cm/s, respectively, leading to a Rossby 
number, R — U/fL , of the order of 10~ 2 to 10“ 3 in mid- 
latitudes. The small Rossby number limit is by itself in- 
sufficient to justify a linear physics on large scales; a ref- 
eree has, in particular, noted that geostrophic motion of 
the second kind, i.e., obeying the classical thermocline 
equations, is intrinsically nonlinear. However, over the 
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Figure la. Domain of the Pacific Ocean model and 
the number of model layers which reflect the topogra- 
phy used during the estimation studies* Vertical extent 
of layers is 100, 500, 1000, and 2400 m (from top to 
bottom), giving a maximum model depth of 4000 m in 
the center of the figure. 

timescales we are treating, the oceanic response appears 
to be dominated by near-adiabatic, wave-like motions 
which are linear. Our fundamental justification is that 
the method being used appears to “work” in the sense 
that the principle of superposition applies for scales and 
amplitudes of the observed differences and for periods 
up to a year. 

We will postulate, in addition, that modeling skill 
will continue to evolve so that even if some present gen- 
eration model fails the test of closeness to the data, 
future models will eventually pass it. (By this means, 
we achieve a separation of the most immediate prob- 
lems involved in model construction from those arising 
in devising estimation procedures.) 

The model will be forced to consistency with observa- 
tions on large scales, and there is an additional assump- 
tion that such a model will then also have greater skill in 
its computation of the smaller scales (e.g., mesoscales) 
through the scale coupling present in any model, even 
a linear one. Alternatively, one could force the model 
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Figure lb. Sea level elevation (in centimeters) and 
surface current fields of the model after a 23-year spin- 
up phase. Every second velocity vector is plotted. Bold 
and thin lines indicate positive and negative sea level 
elevation with a contour increment of 10 cm. 


to consistency with the small scales and anticipate im- 
proved skill on the large scale. The actual extent 
to which OGCMs are controllable through boundary 
observations confined to restricted ranges of scales is 
largely unknown. One seeks ultimately to use the al- 
timetry at all scales available from alongtrack data; we 
confine the present discussion to long wavelengths be- 
cause it provides an adequate basis to test and demon- 
strate the method without straining the available com- 
puter resources. (There is a resemblance in our strategy 
to that used by Fukumori and Malanotte-Rizzoli [1995] 
in a different dynamical setting.) 

The paper is organized as follows: The model setup 
for the Pacific basin is described in section 2. Section 3 
discusses the methodology of estimation using altime- 
ter data with pressure Green’s functions. Their initial- 
ization and successive evolution are discussed in some 
detail in section 4. The state estimation method is il- 
lustrated in section 5 using synthetic data to provide 
insight into the potential of reconstructing the oceanic 
state (three-dimensional flow field and tracer distribu- 
tion) from surface height/pressure observations. Real 
T/P sea surface height observations are used in section 
6 to infer the present-day Pacific circulation and mass 
distribution. 

2. Model Configuration 

The model used in this study is the primitive equa- 
tion Geophysical Fluid Dynamics Laboratory (GFDL) 
model [Bryan, 1969; Cox , 1984], implemented for 
present purposes in the Pacific Ocean north of 30°S 
with realistic coast lines and bottom topography (Fig- 
ure la). Because the GFDL model is widely used, only 
those elements necessary to understand the present con- 
figuration will be summarized here. 

The underlying governing equations which make use 
of the Boussinesq, hydrostatic, and traditional approxi- 
mations are conservation equations for momentum, vol- 
ume, potential temperature 9 , and salinity 5, aug- 
mented by a diagnostic equation for density. Bottom 
and side walls are insulating, i.e., [dn(0, S) = 0]. A 
no-slip side wall condition is imposed and the bot- 
tom is “free-slip.” At the surface there is a rigid lid, 
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Figure 1c. Surface temperature field corresponding to 
Figure lb. 
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Figure ld,e. (d) Coarse 10° by 10° grid, on which the 
Green’s functions where calculated in layers 1 through 
3. (e) Owing to geographical shoaling, the bottom layer 
has a smaller geographical extent on the coarse grid. 

and the momentum flux is given by the wind stress 
[d z (u, u) a with A, <f> being longitude and lat- 

itude, respectively. Surface heat and freshwater (salt) 
fluxes are mimicked by Newtonian damping terms, with 
the model 9 and 5 fields relaxed toward Levitus’s an- 
nual mean surface ^L ev , Sl cv fields as H e = 7(^Lev 
and Hs — T^Lev “ $) with timescale 7 ' 1 . 

Because our main intention here is the exploration of 
the method in the context of real data, a limited model 
resolution is adopted, with a 1° horizontal grid spac- 
ing and with only four layers in the vertical of thickness 
100, 500, 1000, and 2400 m (from top to bottom). The 
maximum model depth is 4000 m. The southern wall 
and the Indonesian passages are artificially closed for 
this specific realization. Although this configuration is 
somewhat limiting in providing a highly realistic model 
of the circulation of the Pacific, it significantly simpli- 
fies the situation without changing the fundamental na- 


ture of the problem. In future applications, the resolu- 
tion will be increased substantially, and ultimately, the 
model will have a global domain. 

Starting from Levitus’s [1982] annual mean 6 and S 
distributions and a resting flow field, the model was in- 
tegrated forward for about 18 years (160,000 time steps) 
with a mean surface wind stress provided by Trenberth 
et a L [1989] from routine analysis of the European Cen- 
tre for Medium-Range Weather Forecasts (ECMWF) on 
a 2.5° x 2.5° grid and averaged over the period 1980 to 
1986. The time step At was 1 hour, and parameters 
for mixing and diffusion were used as listed in Table 1. 
Owing to the relatively strong surface T and S forcing 
and its damping effects on the surface pressure Green’s 
functions (see below), the relaxation time scale 7 " 1 was 
increased from its initial 30 days to 100 days, and the 
model was integrated forward another 4.5 years (40,000 
time steps), giving a total of roughly 23 years of spin-up. 

The surface pressure and velocity fields at the end 
of the 23- year spin-up are shown in Figure lb. The re- 
lated surface temperature field is shown in Figure lc. In 
the North Pacific, the major circulation components are 
simulated. However, because of the low vertical resolu- 
tion, there is not much resemblance to observations in 
the tropics, and the circulation south of the equator is 
dominated by the artificially closed southern boundary. 

3. Methodology 

The estimation procedure introduced below is for- 
mulated for observed model/data differences in surface 
pressure. Our goal is to describe those differences as 
a solution to linear model dynamics in terms of model 
Green’s functions, using all observations available dur- 
ing a period r. As explained below in detail, the result 
is a linear inverse problem for a set of disturbance coeffi- 
cients, which is then solved. Generally, the linear model 
can be set up on the full OGCM grid. (To avoid ambi- 
guity, the original model will henceforth be referred to 
as the “OGCM.”) However, consistent with the scale ar- 
gument given in the introduction, the linear estimation 
problem is defined on a coarse grid, obtained by dividing 
the ocean into a series of three-dimensional rectangular 
boxes of horizontal scale L — 10° with spatial cover- 
age as shown in Figure Id. A stack of such coarse grid 
boxes occupies the whole water column with vertical 


Table 1. Model Parameters 


Parameter 

Symbol 

Value 

Resolution 

Horizontal 

A<t> } AA 

1° 

Vertical 

A 2 

100, 500, 1000, 2400 m 

Maximum depth 

Ho 

4000 m 

Time step 

At 

1 hour 

Relaxation coefficient 

A- 1 

30/100 days 

Horizontal mixing 

•^hm 

10 4 m 2 s' 1 

Vertical mixing 

■Avm 

10~ 2 m 2 s _1 

Horizontal diffusion 

^hh 

10 7 m 2 s _1 

Vertical diffusion 

^vh 

0.5 x 10~ 4 m 2 s~ l 
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dimensions defined by the layer thickness. Boxes next 
to the lateral boundaries are not included in the set of 
perturbed regions because the presence of intense cur- 
rents and other boundary effects can more readily lead 
to violations of linearity. The bottom layer (Figure le) 
has a smaller geographical extent owing to topographic 
shoaling. Only disturbances larger than a horizontal 
scale L will be examined in this paper. However, in 
more elaborate applications, this limitation will be re- 
laxed, permitting smaller scales in the estimation pro- 
cess. Ultimately, one needs to take all scales into con- 
sideration. The result is a linear model, which is much 
reduced in dimension from that of the OGCM, both be- 
cause of the coarser spatial resolution of the former, but 
also because the linear model physics are now embodied 
mainly in the (stored) Green’s functions rather than in 
the full model state vector itself. Achieving this state 
vector dimension reduction is one of the chief goals of 
the present method. 

For notational simplicity, we will apply a one-di- 
mensional counting index to the full three-dimensional 
grid. In addition, t = {0, 1, 2, • * * , tj} denotes a dis- 
crete time index with t j defining the period of inversion 
tj — tj • At, and At being the time step of the linear 
model. At can differ from the time step of th£ OGCM 
and is chosen here to equal the 10-day sampling period 
of T/P. The final tim etj corresponds to a 1-year period 
of T/P data. The sensitivity of the results to the linear 
model time step and the length of the data period are 
subject to further investigation. 

To establish notation, we write the linear homoge- 
neous model in canonical form as 

P {t + 1) = r{t)p{t) (4) 

with initial conditions p(t = 0) = po- Here Y(t) is a 

matrix operator, and p(£) is the oceanic state vector at 
discrete time t = {0, 1, 2, • • •}. Given a pressure Green’s 
function, G p (t,£'), defined from 

G*(t+M') = T(*)G*(t,*'), ( 5 ) 

with initial conditions G = 6 ti t f {6t,t' being the 

Kronecker symbol), a solution to (4) can be written as 
[e.g., Roach , 1982] 

P(*) = I>o(t') GP (M')- (6) 

V 

We use the superscript p, denoting the pressure Green’s 
function, because later we use the associated tempera- 
ture, salinity, and velocity disturbance fields, to be de- 
noted G 0 , G s , and G v , respectively. 

For the state estimation procedure we make use of 
(4)-(6). For this purpose, a linearized version of the 
fully nonlinear model is derived empirically, by disturb- 
ing the spun-up, and forced OGCM by a small pertur- 
bation. The subsequent model evolution (over a 1-year 
period) relative to a similar, but unperturbed, parallel 
run defines the linear model according to (4), with the 
state vector p(i) denoting the difference in the states 


of the two parallel OGCM runs. Accordingly, we de- 
rive the impulse response function, G p { - (£, t'), of the 
linearized model as the time-evolving response of the 
perturbed OGCM relative to an unperturbed parallel 
run, at any grid point i of the full three-dimensional 
OGCM grid and at discrete time t, with a unit ampli- 
tude pressure anomaly imposed at grid point j and at 
t — t'. By definition, G p - (t\ t f ) = 1, and 

= (7) 

the latter condition stating causality, i.e., the model 
will not have responded to future perturbations. Here 
and below, the indices i denote the grid points where 
a response is observed, while the j grid corresponds to 
points at which perturbations are permitted. (There 
is some ambiguity introduced by our use of the label 
“Green’s functions.” We do not expect our empirical 
linear model to be self-adjoint, and for such systems 
a second set of, adjoint, Green’s functions can also be 
defined [Morse and Feshbach , 1953]. Our definition of 
G\ j(t t i') is that of a unit impulse response function for 
the linear model itself as defined in (6).) 

Note that the linear model is obtained here through a 
“black-box” procedure, in which the impulse responses 
are found empirically, assuming that the OGCM will 
respond in a linear fashion to small, local perturba- 
tions. However, we also will make use of the simplifying 
assumption that the Green’s functions are translation 
invariant in time, that is, G p (t,t / ) = G p (< — £',()). If 
the model basic state changes significantly, e.g., if there 
are large seasonal fluctuations, this assumption will fail, 
and we would be forced to compute different Green’s 
functions for different seasons. The major drawback 
to such a need would be greatly increased storage re- 
quirements. Within the limits of the available altimeter 
data, we have found no indication that this refinement 
is necessary. 

We assume now that Green’s functions (£,£')} 
have been determined for a complete set of j points, 
j =r {1,2,- •*,-/}, and shift to the estimation problem. 
In the present context of altimetric sea surface height 
data, observations are restricted to grid points residing 
in the coarse resolution surface layer S C 5, where 5 
denotes grid points in the surface layer and B is the 
three-dimensional grid. Then 

&Pi{t) = (ffPCt(0 - Pi(0) * 6 5 (8) 

defines the difference in surface pressure between the 
OGCM first estimate pi(t) and the observed value 
gp(i(t). The first estimate Pi(t) is obtained by inte- 
grating the OGCM forward in time starting from the 
given initial state with prescribed boundary conditions. 

In the absence of any of the errors discussed below, 
the f*Pi{t) defines the “observations” for our linearized 
estimation problem. However, in practice, there are al- 
ways errors present, and instead the observations must 
be written as 

yi(t) =6pi(t) + ni(t). (9) 
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According to the definition of Sp^t), the term rii(t) 
arises from errors in both the data and the model (see 
below). We seek to find the best estimates a of the 
true forcing coefficients a and their corresponding un- 
certainties P a ,a ( see (20)), which lead to a best fit <5p 
of the true 6p in an optimal sense, 

j = l t '= 0 
J b 

= 0 ). ( 10 ) 

j = 1 t'=0 

plus an estimate of the yet to be determined residual 
(noise) term hi(t). This state estimation process can 
be interpreted as the projection of the noisy “observa- 
tions” onto the dynamical set of linear basis functions; 
the Green’s functions. In (10), the second equality fol- 
lows from translation invariance in time. At the final 
time step t/, the convolution (10) is carried over the 
entire period for which data are available. However, if 
a disturbance in the ocean were fully dissipated after, 
e.g., 50 days, then tj could be reduced, accordingly. In 
our application, tj is assumed to be the whole 1-year 
period of T/P data. 

An estimate Pi(t) of the true absolute pressure Pi(t) 
at any location of the full model domain i and at time 
step t } consistent with altimetric surface observations, 
is then obtained from the sum 

Pi(t) = Pi{t) +<5p.W- (11) 

We pause here to introduce an example: suppose, hy- 
pothetically, the OGCM predicts the correct pressure 
at t = 0 and subsequently everywhere in the model do- 
main, except for the contribution owing to an unknown 
perturbation at grid point jo at £ = 0, which subse- 
quently propagates through the model area. Such a 
perturbation could, e.g., be associated with an error in 
the wind stress field, used to derive the OGCM. Then 
determination of a jo (0) properly initializes the model, 
while all a j0 (£'), t' > 0, and all a ; (£'), j ^ jo, at all t ! , 
vanish. In this special case, (10) would read 

6p l (0=a JO (0)Gf Jo (t ) 0), (12) 

which provides a perfect estimate of the anomaly every- 
where for all future times, and the sum (11) would re- 
produce the observations with no further correction re- 
quired. A further unpredicted disturbance at the same 
grid point jo and t — 1 would require the determination 
of the two unknowns a ;o (0), d ;o (l), leading to 

6Pi(t) = &io(0 )G p ijo (t, o) + a> 0 (l)G? #0 (t, 1). (13) 

In general, one solves simultaneously for all d ; (t') using 
all the observations, and summing over all j produces 
the general linearized response (10) at grid point i £ B. 
The collection {c*j(t')} represents the state vector of 


the linear model. The assumption of linear dynamics 
made in (8)-(10) is that disturbances are superposable 
to adequate accuracy. 

The process of model correction in (10) can be 
thought of as one of constantly reinitializing the lin- 
ear model to the extent that unpredictable forces act, 
or disturbances propagate into the estimation domain. 
The determination of the required perturbation values 
(corrections) is made in such a manner that the result- 
ing model trajectory best fits the observations over the 
entire data period and model domain in a way consis- 
tent with the linear dynamics inherent in the Green’s 
functions. This procedure contrasts with the so-called 
“blending technique” [e.g., Ghil and M alanotte- Rizzoli , 
1991], by which models are re-initialized by direct in- 
sertion of the most recent observations. 

In our context, T/P observations are prescribed only 
on the coarse grid in the surface layer i = {1, 

To render the estimation problem at its simplest level, 
j could similarly be confined to the coarse grid sur- 
face layer, i.e., the special case in which external dis- 
turbances are confined to the surface. In either case, 
the collection of observations over the complete data 
time span form a set of simultaneous equations for the 
unknown coefficients a 


- 0^(1, 0) ... g? i7 (i,0) 
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... 0 


G5.i(1.0) ••• GJ,(1,0) 
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yi(b) 



Q j(t, -1) 


. Mb) . 


_ 

yi(b) . 

(14) 

In matrix form, (14) 

is 







G p a-fn = y 




(15) 


where n represents the noise in all positions at all times 
and y respresents the collection of observations. How- 
ever, more generally, one must accommodate distur- 
bances at depth, which can arise both from initialization 
errors and from errors propagating into the estimation 
domain. Then (10) is generalized to having the j grid 
distributed over the full three-dimensional domain with 
the surface observations being a summation of contribu- 
tions from perturbations occurring in all model layers. 
The Green’s function is then the response ob- 

served in the surface layer at grid point i due to a pres- 
sure disturbance imposed at time step t f at any surface 
or subsurface location. 

In an equivalent formulation, which provides further 
insight into the vertical structures giving rise to ob- 
served surface perturbations, one would rewrite each 
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Green’s function as a set of vertical basis functions. 
Suppose, for example, that the ocean has a flat bottom. 
Then the linear dynamical modes m = {0,1,2 , ..M} 
would provide a natural choice of a vertical basis. Dis- 
turbances at point i would then be composed of pertur- 
bations by mode m ) and resulting anomalies 6p m (t) 
would be the history of the pressure anomalies of each 
mode m. Empirical orthogonal functions (EOFs) or 
any other basis could be used similarly. However, be- 
cause linear dynamical modes in the presence of varying 
bathymetry are defined only locally, we use formulation 
(15), in which each of the coefficients excites energy in 
the vertical modes. Including the subsurface coefficients 
permits one to construct the equivalent of modes, either 
barotropic or baroclinic. 

The system is a linear time-evolving one, and there 
are several options for determining the state vector 6l. 
To render the discussion simple, we will consider a solu- 
tion obtained by the method of singular value decompo- 
sition (SVD [e.g., Wunsch , 1996]) which uses the entire 
domain of observations all at once. If (15) becomes too 
large, sequential methods, such as filters and smoothers, 
or recursive methods can be employed. Application of 
these standard approaches would normally lead one to 
use the canonical form (4) ( Y is readily determined 
from G p by (5) [see Menemenlis and Wunsch , 1996]). 

We solve (15) after column scaling, written as 

G p W _ ^W^q + n = y, (16) 

where W has the diagonal form 

W = diag {(• * • h\/ Ho •••),(••• ^ 2/ #0 * * *)> 

(-W- ■■).(• -h 4 /H 0 •••)}> ( 17 ) 

with hi denoting the layer thickness of the ith layer 
and H 0 is the total depth. Let the SVD of the Green’s 
function matrix G P W ~2 be 

G p W" 5=UAV t , (18) 

with U and V containing the left and right singular 
vectors, A having the singular values on its diagonal, 
and the superscript T denoting the transpose. Then 
the estimated state vector a is 

« = W-iVjcAi l U£y, (19) 

with an estimated uncertainty 
P Qa =< (a - a )(a - a ) T > 

= a 3 W-S (VjcA^V^+Q.RaaQ') W"L 

(20) 

distinguishing the true ol from the estimated one d . 
Here K is the rank of the system and is determined 
in practice by prior assumptions about the white noise 
variance a 2 and the magnitude of the singular values 
A i (not to be confused with longitude A); the noise is 
discussed further below. Column vectors of are the 
null space vectors, R aa is the second moment matrix 
of the null space coefficients, and P aa is the sum of 


the noise uncertainty and the uncertainty from the null 
space variance. Information about R aa must be pro- 
vided from prior information; in its absence, one can 
speak only of the unresolved components. See Wunsch 
[1996] for further details. 

Given estimated coefficients a, the linear model vari- 
ables are obtained on the full three-dimensional grid, 
i G B, from equations like (10) as a function of time, 
with 

J h 

6Pi(t) = £ E - t\ 0) (21) 

j V — 0 
J if 

[SOi'SSiMm = E E «i(OGS S,V) (t - 0) (22) 

j t'=o 

where the matrices G^ ,5,v l describe the response of the 
temperature, salinity, and velocity fields to the associ- 
ated pressure perturbation. The velocity field also fol- 
lows readily from the resulting pressure field by geostro- 
phy, and in practice, 6\i(t) was found this way. 

Uncertainties of the estimated fields are obtained 
from 

P pp = G p P aa G pT (23) 

P [ »,s.v]= Gl*- s > v ]p aa Gl e ' s ' v ] T (24) 

A few additional comments seem to be useful. First, 
the linear model trajectory is taken about the prior es- 
timate OGCM state in an assumption essentially the 
same as used in what is known as the “linearized 
Kalman filter” [see Wunsch, 1996]. A generalization 
would be to linearize about a new state (the OGCM 
plus linear model state) which in the sequential es- 
timation context, leads to the “extended Kalman fil- 
ter” and associated smoothers. There is a parallel set 
of linearizations used in adjoint/Pontryagin principle 
methods in which the adjoint model represents a lin- 
earization about the instantaneous OGCM state, which 
will change through the iterative solution used in that 
method. There is also a close connection of the present 
method with that of “representers” (see Bennett [1992] 
for details). As with any filter-smoother approach, (14) 
uses observations from the entire estimation period to 
determine the estimated state at any specific time. This 
can be understood by recognizing, e.g., that a(0) is mul- 
tiplied by G(N , 0) to couple future observations with 
estimates at prior times. 

No attempt has been made here to evaluate directly 
the relative numerical efficiencies of the present “whole- 
domain” approach relative to that of a sequential or 
other solution method (although, in general, whole 
domain methods, because they do not require multi- 
ple computations of the error covariance as filters or 
smoothers do, are far more efficient). Our focus is on 
demonstrating the concept of numerical linearization as 
a workable methodology. In general terms, however, the 
Green’s function approach can be thought of as being 
related to an inverse operator which is vastly reduced in 
dimensions to contain only that kernel required by the 
specific data type and distribution: the surface pressure 
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boundary condition. In contrast, the adjoint operator 
of the (linearized) model would need to invert the com- 
plete model state in space and time. 

It is possible that when Green’s functions solu- 
tions are computed with high resolution, numerical 
pathologies will arise requiring spatial and/or tempo- 
ral smoothing, but we do not yet have experience with 
this situation. 

Because observations are taken here to be model/data 
differences, the “noise” n*(£) represents all the errors of 
the perturbation pressures at the surface, including re- 
maining environmental or geophysical data errors from 
such sources as tides and orbits, and errors due to miss- 
ing model physics and insufficient forcing fields and for 
processes on scales smaller than the dimension of the 
Green’s function representation areas. We will gener- 
ally assume that 

< n(t) >= 0, < n(t)n T (t) >= R (t), 

< n(t)n T (t') >= 0, t / t'. (25) 

The determination of the actual noise level and struc- 
ture must be addressed using real data. Systematic 
model errors will not be adequately represented through 
the joint covariance R, nor will systematic (geograph- 
ically correlated, time-independent) errors in the T/P 
observations. Careful posterior tests are required to test 
the assumptions. Uncertainties in G p are also present 
owing to the initialization procedure and the subsam- 
pling scheme. Discussion of the influence of such errors 
on the results is a nonlinear problem (sometimes called 
“total leastsquares” [e.g., van Huff el and Vandewalle , 
1991; Wunsch , 1996]). We will ignore such errors in 
this present discussion. 

4. Finding Green’s Functions 

Conceptually, we seek to determine the response of 
the spun-up OGCM to an initial, regionalized, and 
weak, surface elevation/pressure perturbation 

C(0? — (o^r,r 0 £t,t 0 (26) 

confined to a single surface grid point at the geographi- 
cal location ro = (<f> o, A 0 ) in the surface layer. However, 
with a rigid lid, the surface elevation is not part of the 
OGCM state vector, and we cannot directly perturb the 
model with a surface pressure perturbation. To within 
the quasi-geostrophic approximation, however, the re- 
sponse to an initial geostrophically balanced vortex, 

v = y(k x V()e"(^^) , (27) 

again confined to a single surface point as £(<£,A) = 
Co exp {- [(A - \ 0 ) 2 + (<j> - <f> 0 ) 2 ] /L 2 } (Figure 2), is 
equivalent to imposing a pressure perturbation on the 
model. The spacescales and timescales are extended 
to L = 5° and r — 5 days to represent an initial 
pulse-like perturbation after its geostrophic adjustment 
process on a rotating sphere. The general problem of 
geostrophic adjustment has a large literature going back 



Figure 2. Schematic of the model perturbation. 
An initial height disturbance is transformed into a 
geostrophically balanced vortex, which is imposed on 
the model surface layer. See text for details on 
spacescales and timescales of the initialization proce- 
dure. 


at least to Hough [1897], and including the famous pa- 
per of Rossby [1936]. Blumen [1972] and Gill [1982] 
provide comprehensive reviews. From those studies, an 
initial adjustment timescale of the order of several in* 
ertial periods is found. 

Green’s functions were computed and stored for ini- 
tial disturbance (27) in each of the coarse grid points. 
The linear model pressure Green’s function G p is ob- 
tained by the normalization of the resulting time series 
of the pressure perturbation by its maximum value (re- 
sulting in G v - (0) — 1). The matrices G^ ,s,v l are ob- 
tained from the related temperature, salinity, and ve- 
locity fields. 

The adjustment process of the interior ocean to sur- 
face pressure disturbances is of interest in its own right. 
(One must distinguish perturbations within the upper- 
most layer of the ocean from those imposed as atmo- 
spheric loads; for the latter, see Ponte ) [1992], and Wun- 
sch and Stammer , [1996].) We confine the discussion 
here to a brief summary of the response as obtained 
from the numerical model. Analytical solutions of rele- 
vance to our computations can be found in the works by 
Bolin [1953], Fjeldstad [1957], Longuet- Higgins [1965], 
Geisler [1970], Philander [1978] and Blumen [1972]. In 
a different context, a related model response obtained 
from a subsurface density perturbation was analyzed by 
Sarmiento and Bryan [1982], which in its mass trans- 
port stream function shows a wavelike pattern in the 
central North Atlantic, very similar to the results dis- 
cussed below. 

The response of the model to a localized surface pres- 
sure disturbance associated with the isolated vortex 
(27) centered at 35°N, 220°E is depicted in Figure 3 in 
terms of its surface pressure and flow field after t = 16 
and t — 41 days (from now on, we use t as dimensional 
time). An expected asymmetric pattern of a radiating 
Rossby wave field is visible, which rapidly spreads over 
the entire basin and reaches the Kuroshio within about 
10 days. The dominant wavelength decreases with time, 
as the longer waves are reflected and dissipated. Af- 
ter about 50 days, the barotropic wave field is largely 
gone. The influence of topographic features such as 
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Figure 3. A Green’s function, illustrated from the surface pressure and velocity fields (top and 
bottom panels, respectively) at t — 16 (left) and t — 41 days (right). The initial perturbation 
was imposed at 30°N, 215°E. Note the pronounced influence of the bathymetric features on the 
propagation characteristics of the developing fast Rossby wave field. Maximum values are clipped 
in Figure 3a. The contour increment is 0.1 cm. 



Figure 4. Same as Figure 3, but for a perturbation imposed at 5°N, 215° W and subsampled 
at t - 25 days (left) and t = 58 days (right). Unlike the perturbation in midlatitudes, the 
enhancement of a low-mode Kelvin wave at the equator and its reflection at the eastern side is 
visible. The contour increment is 0.1 cm. 
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the Hawaiian Ridge or the Emperor Seamount Chain 
is apparent in Figure 3b, effectively blocking the fast 
barotropic waves. 

The model response, which in midlatitudes shows 
predominantly westward propagation of locally imposed 
signals, is altered in the tropics: a velocity anomaly im- 
posed at 5°N, 215°E (Figure 4) generates a fast east- 
ward moving low-mode baroclinic Kelvin wave along 
with a westward going baroclinic equatorial Rossby 
wave. The Kelvin wave is partially scattered by the East 
Pacific Rise before reaching the eastern boundary. In 
agreement with the theory of equatorial dynamics [e.g. , 
Philander , 1978], the low-latitude disturbance, being 
predominantly baroclinic in vertical structure, spreads 
zonally but does not penetrate far into high latitudes. 

In the vertical, the model response is generally a sum 
of all dynamical modes (the barotropic plus three baro- 
clinic) which are excited by the initial disturbance. Fig- 
ure 5 depicts a zonal section of the meridional velocity 
component through the center of the initial disturbance 
over a 1-year period. As can be seen from the figure, 
the fast initial barotropic response is followed by the 
emergence of increasingly high baroclinic modes, each 
of which moves westward with its associated Rossby 
wave phase speed. The anticyclonic vortex imposed ini- 
tially only in the surface layer induces an ageostrophic 
flow component into the vortex center and an associated 
secondary vertical velocity field (not shown) with down- 
welling in the vortex interior and upwelling at the outer 
edges which is associated with the return flow in the 
vertical. The consequent shift of the isotherms leads to 
a warming in the vortex center and a cooling at its outer 
edges. The complex interior response is important in a 
general sense: surface elevation/pressure is coupled to 
interior motions over the entire water column, motions 
which we can hope to infer from an inverse computa- 
tion. 

When a perturbation is introduced in a subsurface 
layer, a number of expected changes can be observed. 
Figure 6 shows the response to a perturbation in either 
of the four layers at 35° N, 205° E after 83 days. The re- 
mote barotropic field is largely independent of the spe- 
cific depth of the initial disturbance. In the nearfield, 
the details of the baroclinic response are different. Con- 
sistent with the dominance of the first baroclinic mode, 
the immediately overlying surface response to a pres- 
sure disturbance in layers 3 and 4 is of opposite sign, in 
contrast to a perturbation imposed in the two surface 
layers. In addition, details of the energy partitioning 
of baroclinic modes depends on the depth of the initial 
perturbation (e.g., giving rise to enhanced second-mode 
energy in low latitudes when initialized at middepth). 

5. A Twin Experiment 

We now turn to the estimation problem: to calculate 
the full three-dimensional oceanic state at all times of 
interest as outlined in section 3, given only altimetric 
observations of the sea surface elevation field (or equiv- 
alently, the related geostrophic flow in the surface layer) 
and a model forecast. 


Although our primary interest is in the results from 
actual data, we nonetheless will digress briefly to de- 
scribe the results obtained with perfect “data” from a 
so-called twin experiment. Real observations raise ques- 
tions which are rarely encountered in discussions of es- 
timation with artificial data, but use of artificial data 
permits us to isolate problems with the method which 
are independent of the complexities of real observation 
uncertainties. The purpose of this intermediate step is 
to establish the estimation problem in the presence of 
a completely known ocean state. Because of the ab- 
sence of good in situ data coverage simultaneous with 
the altimetry, a twin experiment is the only basis for 
exploring a comparison over the full depth range. 

Starting from the spun-up state, the OGCM was inte- 
grated forward for one additional year, the time history 


122* E 142" 162® 1B2® 202® 222® 242® 262® 202° E 



Figure 5. Zonal sections of the northward (u) veloc- 
ity component (corresponding to Figure 3) through the 
center of the initial velocity perturbation for (a) t — 16, 
(b) 41, and (c) 410 days. Negative values are stippled. 
An arbitrary contour increment is shown, to indicate 
vertical structures. 
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Figure 6. Sea surface pressure fields, which result from putting velocity perturbations into layer 
(a) 1, (b) 2, (c) 3 and (d) 4, after 83 days. The contour increment is 0.05 cm. 


of which gave the first-estimate reference state. To gen- 
erate artificial “model data,” the OGCM forcing was 
changed from constant to monthly mean wind stress 
values, and run for an additional 10 years into a new cli- 
matological equilibrium. The subsequent eleventh year 
was taken as the “true” state, which is supposed to be 
recovered from the surface “observations,” in our frame- 
work, taken to be the surface elevation anomalies refer- 
enced to the previous reference state. As an example, 
Figure 7 shows instantaneous anomaly fields of surface 
pressure and velocity, representing winter conditions. 
Such anomaly fields were sampled on the coarse grid 
every 10 days over a 1-year period to form the vector of 
observations y as stated in (14). 

To obtain more insight into the problem of estimating 
the full ocean state from surface height/ velocity obser- 
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Figure 7. Instantaneous difference fields of p s and hor- 
izontal surface velocity v between the synthetic “model 
data” and the OGCM first guess, representing winter 
conditions. The contour increment is 10 cm. 


vations, various experiments were performed, all listed 
in Table 2. In a first experiment, labeled Gl, both ob- 
servations and coefficients a are confined to the coarse 
grid in the surface layer only. There are then 2590 un- 
known values (74 points times 35 time steps of 10 days 
each), and an equal number of observations connected 
by (14). Confining the unknowns to the surface layer 
is equivalent to the assumption that only in the surface 
layer do unpredictable changes take place, e.g., because 
the wind perturbed the system. Changes in the thermo- 
cline are then uniquely coupled to surface expressions 
through vertical dynamical modes. Knowledge of the 
forcing coefficients in the surface layer are assumed to 
completely control the deep ocean. 

In a second experiment, labeled G2, the a are allowed 
to be nonzero in all four layers over the full data period, 
thus producing 8960 unknowns in 2590 equations (15). 
Here we make the extremely pessimistic assumption 
that disturbances leading to later changes in surface 
pressure can appear at depth within the model domain 
from unknown causes, e.g., missing model physics. It is 
plausible that such perturbations can enter the domain 
at depth across the lateral edges, but permitting them 
also within the interior represents an extreme case. 

Two further experiments were performed which are a 
combination of the previous two scenarios. The first of 

Table 2. Definitions of the Four Twin Experiments 
Exp. Comment 

Gl a / 0 only in surface layer for all t 

G2 a ^ 0 in all four layers for all £ 

G3 same as Gl, plus a ^ 0 in all layers at t = 0 

G4 same as Gl, plus a ^ 0 in all layers at t = 0 and 

along lateral boundaries for all t 
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Figure 8. Singular value distribution from the two 
cases G1 and G3. See text and Table 2 for details on 
the experiments. 

these, G3, is basically experiment Gl, but with nonzero 
ct(0) in the full water column, equivalent to permitting 
an initial model correction over the full water column. 
Later developments are, however, through surface forc- 
ing perturbations alone. Finally, G4 includes in addi- 
tion to G3, all a(t'), t f > 0, along the lateral boundaries 
of the coarse grid, thus accounting for perturbations 
which might enter into the estimation domain at depth 
from outside the domain. 

The singular values A* of the two cases Gl and G3 
are depicted in Figure 8. Formally, both systems are 
full-rank (the same being true for the other two cases), 
with no singular values actually vanishing. For Gl, the 
full rank means that the solution is completely resolved 
and would be perfect if the synthetic observations were 
actually perfect. There are five very small A*, how- 
ever, and some rank deficiency is expected. Errors do 
arise from the loss of spatial resolution in moving to 
the coarse grid, and from edge effects in the averaging. 
When the five very small A, and corresponding singular 
vector structures are dropped, some data and solution 
resolution is lost, basically confined to the tropics, and 
mainly at the beginning and end of the 1-year time pe- 
riod. A somewhat more realistic situation is obtained 
by taking the effective rank as K = 1500, which im- 
plies a noise level of about 3% of the surface elevation 
variance. (Results are not very sensitive to the specific 
rank. The value of K — 1500 was chosen here for con- 
sistency with the application to T/P data described in 
the next section.) 

For G3, and similarly for G2 and G4, the A; vary only 
over about 2 orders of magnitude. Whether full rank 
or not, cases G2 to G4 are always formally underdeter- 
mined, and one normally expects to use prior statistical 
information and understanding in finding a best solu- 
tion. In the spirit of exploring a somewhat pessimistic 
situation (G2), we assumed only that the solution a 
has a variance inversely proportional to the layer thick- 
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Figure 9a- Diagonal elements of the data resolution 
matrix T u = UU T , which results from a rank reduction 
to k — 1500 and plotted for t — 100 days in geographical 
order. 

nesses (smallest in the deep water). For G3 and G4, we 
altogether suppress any significant variances in in 

the interior, except at t l — 0 for an initial correction 
(G3) and along the lateral edges in layer 2 and 3 for 
t f > 0 for G4. A useful refinement (not done here) would 
be to impose larger variances near western boundaries 
or near mean currents. The imposed variances could 
come from the model itself, or from observations of any 
kind. 

Gl and G3 produce basically identical results, but 
the latter case is more physically attractive and we have 
chosen it as the standard “best case” to be described in 
more detail. Differences in the other experiments and 
the physical implications are discussed subsequently. 

In the present artificial situation, the main issues are 
those of resolution. A typical mid-data-stream spatial 
pattern (corresponding to t — 100 days) of data and so- 
lution resolution is displayed in Figure 9 on the coarse 
grid. The fields are obtained assuming an effective rank 
of K — 1500, in which modes with small spatial scales 
and timescales are mostly suppressed. Generally, the 
data resolution (Figure 9a) is enhanced in higher lat- 
itudes and is small along the equator, with minimum 
values residing along 5°S in the eastern basin. In addi- 
tion, some fine structure is visible, which to a limited 
degree shows a correlation with the bathymetry, e.g., 
along the Emperor Seamount Chain and the Hawaiian 
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Figure 9b. Same as Figure 9a but for solution (param- 
eter) resolution matrix T* = VV T in the first model 
layer. 
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Figure 10. Time means of (left) <5p(t) and (right) 6p 0 bs{t) in centimeters in top three model 
layers (from top to bottom) at all 74 locations on the coarse grid. 


Ridge. The fact that all values are significantly lower 
than unity indicates that no individual datum is fully 
resolved, i.e., the system lacks adequate information to 
distinguish individual equations from linear dependence 
on one another [ Wunsch , 1996]. Owing to the presence 
of fast modes on large zonal scales, in the band ±10° 
around the equator, individual observations sampled at 
a 10 day interval are barely independent of each other. 
In particular, data from the eastern tropical Pacific are 
least important in constraining the solution. One physi- 
cal reason is the presence of the East Pacific Rise, which 
hinders the penetration into the interior basin of the 
equatorial Rossby waves excited in that area. In addi- 
tion, any Kelvin wave generated there soon disappears 
at its eastern border. 

The solution (parameter) resolution for the same time 
step in the surface layer is shown in Figure 9b. Al- 
though slightly different in detail, it shows basically the 
same spatial pattern and amplitudes as appear for the 
data resolution. The similarity suggests a strong, but 
incomplete, dependence of the surface layer values of 
6l on the local measurements. In general, d are deter- 
mined as linear combinations of the elements of the true 


value a, owing to a less than full rank system. Again 
the tropical area shows the largest deficiencies. 

Based on the artificial model observations and the 
SVD of the Green’s function matrix, corresponding co- 
efficients d were obtained from (19) with K = 1500, 
and an estimate of the full model state was recovered 
subsequently from (21) to (22). 

A comparison of the estimated pressure anomaly field 
6 p with the “true” anomalies is shown in Figure 10 for 
the top three layers. Because the bottom layer is only 
partially covered by the coarse grid, it is not displayed 
here; results, however, are qualitatively very similar to 
what is found from layer 3. Shown are instantaneous 
fields from t — 100 days, which nonetheless represent 
typical situations. 

In summary, “true” model fields are simulated in the 
top two layers, qualitatively and quantitatively. The 
rms misfits are reduced to about 10% of their original 
values in the top two layers (Figure 11a), with the misfit 
in layer 2 slightly exceeding that of the top layer. Cor- 
responding cross-correlation coefficients are as high as 
0.95 or above, right from the beginning and over most 
of the 1-year period. Although the cross correlation is 
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Figure 11a. The rms residuals from the twin experi- 
ment at K - 1500 from top three model layers evalu- 
ated on all 74 grid points for each of the 35 time steps 
and normalized by the initial model/ data misfit. Values 
from layers 1, 2, and 3 are drawn by solid, dashed, and 
dotted lines, respectively. 

highly significant in the third layer, growing from an 
initial 0.5 to about 0.8 in the second half of the esti- 
mation period, related rms residuals vastly exceed the 
initial small discrepancy in layer 3. 

A close inspection reveals the troublesome areas 
to be confined to the northern border of the coarse 
grid. There, the Green’s functions, being predomi- 
nantly barotropic in nature, do not project enough en- 
ergy onto the first baroclinic mode, which is necessary 
to match the observations. The tendency is, however, in 
the right direction, with baroclinic modes increasingly 
gaining energy through the continuous reinitialization. 
Improved agreement from assimilation runs lasting for 



Figure lib. Cross-correlation coefficient between 
6p(t) and <5p o bs(0 as a f unc fi° n °f time. Values from 
layers 1, 2, and 3 are drawn by solid, dashed, and dotted 
lines, respectively. 



Figure 12. Space-time mean of the diagonal ele- 
ments of the pressure resolution matrix given as T p = 
G p T t/ (G p ) -1 , averaged over the full coarse grid and 
time span of data for each layer at K — 1500 (dashed) 
and K — 2590 (full-rank, solid). 

longer periods can be anticipated. Furthermore, the 
predominantly barotropic nature of the Green’s func- 
tion in these areas is an artifact of the low vertical reso- 
lution. Another area of potential problems can be found 
in the eastern tropical Pacific, where topographic effects 
are likely to be significant. 

The reduced skill at depth is readily anticipated from 
a resolution analysis of the model and observations, and 
it is the availability of this information which would pre- 
vent any user of the method from falsely inferring that 
the results at depth are better than is warranted, when 
the true field was not actually known. Figure 12 de- 
picts the pressure resolution T p = G P T„(G P ) 1 as a 
function of depth when averaged over the entire time 
span of data and over each layer for K = 1500 and 
K = 2590 (full rank). Here T* is the solution resolu- 
tion matrix as given in Figure 9. As in (21), the ma- 
trix G p is complete, covering the whole space and time 
span, including all four layers, and (G p )“ 1 is the gener- 
alized inverse of this overdetermined system. Not sur- 
prisingly, the solution resolution degrades rapidly with 
depth, even for the full-rank case. Physically, the re- 
sult means that pressure anomalies at great depths have 
unobservable consequences at the sea surface over time 
spans of 1 year. Although longer records will improve 
these results, as deep anomalies generate observable sur- 
face changes over time, it was the anticipation of this 
situation that led Munk and Wunsch [1982] to propose 
the complementarity of altimetry and ocean acoustic to- 
mography, the latter providing the resolution at depth, 
which is difficult for altimetry. 

Temperature anomalies 60 estimated from (22) are 
compared with the “true” fields in Figure 13, again 
taken from t— 100 days and representing typical instan- 
taneous situations. The poorest visual agreement be- 
tween the inferred and true temperature fields is found 
in the surface layer, where the true model state is re- 
laxed continuously toward Levitus’s climatology. The 
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Figure 13 . Time mean (left) 6T(t) and (right) ST Q ^ s (t) in degrees Celsius in top three model 
layers (from top to bottom) at all 74 locations on the coarse grid. 


rms difference is about 1°C (Figure 14a), and the cross 
correlation is as low as 0.4, yet significant, over the en- 
tire period (Figure 14b). We infer that the poor agree- 
ment indicates the inconsistency of the Levitus numbers 
with the model physics in the surface layer. In con- 
trast, general agreement of the results with the “truth” 
is found in the second layer, where the mean correla- 
tion is 0.6. As in the top layer, the largest residuals are 
confined to the tropical ocean, where the induced un- 
certainties are largest. As a result, the temperature is 
systematically underestimated by 2° to 3°C (rms mis- 
fit is again close to 1°C, globally). In the third layer, 
the correlation between observations and estimations is 
only marginally significant. Again an estimate with in- 
creased accuracy in the deep ocean can be anticipated 
from longer data sets because the dynamical adjustment 
timescale in the deep ocean is 1 to 2 orders of magnitude 
larger than near the surface. 

Summarizing the results from the various case stud- 
ies, no significant difference in the estimated state mea- 
sured in rms residual amplitudes of the pressure and 
temperature fields was obtained from G1 compared to 
G3. The effect of permitting c*(0) ^ 0 in all layers leads 
only to slightly decreased residuals in the top two layers 


during the first 50 days, but a small degradation below. 
Permitting subsurface perturbations over the full time 
span during G2 reduces the residuals in surface pres- 
sure somewhat, consistent with the increased degrees of 
freedom. In the second layer and below, which are not 
constrained during the inversion, however, a significant 
degradation of the estimated state relative to the truth 
was found as compared to the purely surface driven 
case. Including only coefficients along lateral bound- 
aries (G4) did not show any benefit near the surface (as 
compared to G1 and G3), but led to some degradation 
in layers 3 and 4. 

6. TOPEX/POSEIDON Altimetry 

Our central goal is the estimation of the present state 
of the ocean given real altimeter data, and in the fol- 
lowing, we will discuss the results obtained with T/P 
data from the first year of the mission. T/P data 
from the 1-year period December 21, 1992, to Decem- 
ber 3, 1993, corresponding to repeat cycle 10 through 
44, were edited and corrected as described by Stammer 
and Wunsch [1994] and King et al [1994]. Only two 
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Figure 14a. The rms temperature residuals from the 
twin experiment at K — 1500 from top three model 
layers evaluated on all 74 grid points for each of the 35 
time steps. Values from layers 1,2, and 3 are drawn by 
solid, dashed, and dotted lines, respectively. 

significant modifications of the standard merged T/P 
geophysical data records produced by the project [Be- 
rta da, 1994] should be noted here because they influence 
the results and the underlying errors. First, to obtain 
observations of the dynamic sea surface height (SSH), 
we referenced the T/P observations to a hybrid geoid 
(JGM-2/OSU91a) derived by Nerem et a l. [1994], and 
Rapp et al. [1991], where the former is used to spherical 
harmonic degree 70 and the latter beyond that to degree 
360. Second, the tidal correction provided by the T/P 
project was replaced by those estimated by Ma et al 
[1994] from the first year of T/P data (version UT/CSR 
1.4). The resulting mean field of this 1-year period is 



Figure 14b. Cross-correlation coefficient between 
60(t) and 60 o b s W as a function of time. Values from 
layers 1, 2, and 3 are drawn by solid, dashed, and dotted 
lines, respectively. 



140“ 160° 180* 200° 220° 240° 260° 280*E 


Figure 15a. Mean TOPEX/POSEIDON sea surface 
height observation inferred from repeat cycle 10 to 44 
relative to the JGM-2 geoid. The contour increment is 
10 cm. 

shown in Figure 15a. A mean difference between the 
model-predicted surface elevation and that observed by 
T/P is shown in Figure 15b. Apart from model er- 
rors, the dominant errors in Figure 15b should be those 
of the geoid. Although geoid errors are geographically 
variable and strongly correlated, for present purposes 
they will be treated as homogeneous and “white.” 

To form a vector of observations, for each individual 
10-day period, instantaneous fields similar to Figure 15b 
were averaged in 10° x 10° areas which coincide with the 
coarse grid. Resulting fields were subsampled along the 
grid points given in Figure Id. The model climatology 
and Green’s functions were used with the T/P data in 
the same way as employed for the twin experiments. 
It should be emphasized that the resulting model/data 
differences are based on the absolute T/P SSH, not just 
the time-dependent part. 

We confine the discussion to G3, which was iden- 
tified in the previous section as being the physically 
most plausible situation. The ability to make inferences 
about the ocean circulation is directly dependent upon 
the noise level in the observations. The discussion of the 
noise level is slightly subtle here because we are simul- 
taneously discussing both the time mean circulation, 
which is sensitive to geoid errors, and the low-frequency 
variability (dominated by the annual cycle), which is to 
first order, independent of such errors. A number of 
studies of geoid error [Nerem et al ., 1994; Tsaoussi and 
Koblinsky , 1994] would lead to a rank K = 500 corre- 
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Figure 15b. The difference of this mean sea surface 
height field relative to the reference model state. The 
contour increment is 10 cm. 
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Figure 16. Forcing coefficients oc(i) in centimeters 
of equivalent surface elevation at various representa- 
tive locations. Shaded areas correspond to the uncer- 
tainty of the coefficients as given by (diag{P aa }) 1 / 2 = 
(<r 2 diag{ Vjf AJ^ 2 } ) 1 / 2 . 


sponding to an error at about 20% of the total variance. 
This rank will be used for studying the time mean. In 
examining the annual cycle, however, the rank is taken 
to be K — 1500, corresponding to a noise level a of 
about 2.5 cm (roughly a 5% error level), which is more 
appropriate to the time-dependent measurements, but 
the time means are removed from these results, as they 
are regarded as noise-dominated. An alternative ap- 
proach would be to completely separate the discussion 
of time-mean and time-dependent fields by generating 
a new data set from the time derivatives of the obser- 
vations. The present method is simpler, however. 

In this preliminary attempt, no use was made of 
known correlations in the error fields (resulting, gen- 
erally in a reduction in the skill of the data compared 
to what is actually possible). Use of such correlations 
is a high priority for future work. 

Consider the higher rank system, K = 1500. 

The time histories of OLj (t f ) and their uncertainties, 
±(cr 2 diag{P aa }) 1//2 , are shown in Figure 16 at a few 
representative locations j. Uncertainties are computed 



Figure 17. Time mean coefficients a as a function of 
geographical positions. 


by taking into account only the first term on the right- 
hand side of (20), omitting the unknown contributions 
from the null space. In general, the largest amplitudes 
occur, as expected, where initial model/data differences 
are large at i! — 0, representing the initialization errors 
with amplitudes gradually decreasing subsequently to- 
ward an asymptotic level. 

Systematic differences between the model and the 
data appear as corrections of uniform sign over large 
areas and long times. Figure 17 shows the average a 
over a full year at rank 1500. As in section 3, we can in- 
terpret these mean coefficients as representing the con- 
sequences of a systematic error in the forcing, e.g., a 
wrong wind stress boundary condition, or incomplete 
model physics. The correction keeps the model tracking 
the observations by depressing the subtropical gyre and 
uplifting it at higher and lower latitudes in a manner 
one anticipates from the mean model/data difference 
(Figure 15b). 

Time histories of the surface pressure corrections 
<5pi(i) are shown in Figure 18 for a few representa- 
tive locations 2 , along with the uncertainty estimates 
(diag{P pp }) 1 / 2 of 6p from (23). As required, the es- 
timates follow the observations in a slightly smoothed 
manner and generally within the uncertainty limits. 

In summary, the method produces reasonable results, 
both visually and in terms of the prior statistical esti- 
mates. As always, however, the degree of agreement of 
model and data is sensitive to the choice of rank, i.e., 
the assumed noise level of model and data. Observa- 
tions are reproduced identically in the surface pressure 
at full rank. Assuming an error of 2.5 cm (A=1500), the 
observed variance and its seasonal variation are simu- 
lated (Figure 19a), with major deviations limited to the 
initial 1-month period. Accounting for geoid uncertain- 
ties (Jf=500), though, leads to a variance level slightly 
lower than the observed mean, and the estimates do not 
show a seasonal cycle. 


195E.45N 195E, 35N 




195E.25N 195E.15N 




Figure 18. Estimates 6p(t) (solid line) and observa- 
tions <$P o b s (0 (dashed line) at the same representative 
positions as in Figure 16. Shading represents the un- 
certainties of 6p(t)i obtained from P pp =: G p P aa G pT . 
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Figure 19a. Variance as it is observed (solid line) and 
simulated at rank K ~ 1500 (dashed line) and K — 500 
(dash-dotted line). 


The asymptotic rms residuals at K — 500 indicate 
that prior error assumptions are consistent with the re- 
sults, leading to a reduction in residuals close to 20% of 
the observed variance (Figure 19b). For K — 1500 the 
error reduction is actually larger than anticipated, lead- 
ing to an rms noise level of 3% (instead of the assumed 
5%). 

6.1. Time Mean Corrections 

Here we focus on the time mean circulation from the 
rank 500 estimate. Figure 20 shows the absolute sea 
surface elevation field in the OGCM before and after 
correction by addition of the linear model estimates. 
The observations and inversion have forced a weakened 
subtropical gyre in the OGCM with its center shifted 
southward by about 10° and with its zonal extent con- 
fined to the western part of the basin, approximately 



Figure 19b. The rms residuals obtained at K ~ 1500 
(dashed line) and K — 500 (dash-dotted line), as com- 
pared to initial model/data misfit (solid line). 



140* 100° 180° 800" 220° 240° 260° 280°E 



140° 160° 180° 200° 220° 240° 200* 2B0*E 


Figure 20. Absolute model sea surface elevation field 
(a) before and (b) after correction by the data. The 
contour increment is 5 cm. 

west of the Hawaiian Ridge. As compared to the T/P 
observations given in Figure 15a, many of the spurious 
small-scale anomalies, which appear to be associated 
with geoid uncertainties, have been reduced by forcing 
dynamical consistency on the motions. 

The vertical structure of the absolute pressure in the 
basic state is shown in the left panel of Figure 21 along 
with the associated geostrophic flow. The time-mean 
correction in pressure to the model first estimate and 
the associated flow are shown the right panel of the 
figure. In general terms, the large-scale gyre is reduced 
from the original OGCM estimate by about 1 cm/s to 
2 cm/s in the top two layers. In contrast, the cyclonic 
deep circulation is increased below 1000 m depth. The 
northern part of the correction has a simple, barotropic- 
like vertical structure, but has a more complex structure 
in low latitudes, associated with the eastward flowing 
equatorial counter current and with weak westward flow 
beneath of it. 

The current and pressure anomalies are associated 
with temperature (density) anomalies at all depths. 
Some impression of the corrections required to these 
fields in the OGCM may be seen in Figure 22 along 
two meridional sections at 165°E and 215°E. Consistent 
with the large-scale changes in sea level, the underlying 
thermocline is cooled or warmed in such a manner that 
the vertically integrated steric changes in density follow 
the required changes in sea level. 

6.2. The Seasonal Cycle 

Because of the heavy spatial filtering, the time- 
varying signal remaining in the T/P data is dominated 
by seasonal changes in SSH with some minor ampli- 
tude fluctuations on timescales of 50 days superimposed 
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Figure 21. (left) Absolute pressure in the basic state along with the associated geostrophic 
flow for the top three model layers (from top to bottom). A spatial mean pressure has been 
subtracted from each field, and positive and negative values are drawn by bold, and thin lines, 
respectively, (right) Time mean corrections in absolute pressure and related geostrophic flow 
which are associated with the estimate given in the left panel. The contour increment is 5 cm. 
The reference vectors represent 10 cm/s (left) and 4 cm/s (right). 


(Figure 18). We therefore focus here on the seasonal 
signal. Figure 23 shows estimates of the surface ele- 
vation anomalies relative to the 1-year mean following 
the inversion for the a at K — 1500. The fields are 
averaged over the four seasons of the year and are es- 
sentially the same as in Plate 3 of Stammer and Wunsch 
[1994]. In contrast to their plate, however, the present 
one is consistent both with the data and with the model 
dynamics. Also shown are the current field anomalies 
associated with the annual cycle elevation corrections. 
Most of these flows are very weak, in agreement with the 
theoretical estimates of Gill and Niiler [1973]. The ex- 
ception is in the tropics, where the seasonally reversing 
North Equatorial Current and the Equatorial Counter 
Current of the order of 2 cm/s are visible. (These 
numbers will increase when the model resolution is im- 
proved.) 


The vertical structure of the seasonal current anoma- 
lies can be seen in Figure 24 in a section of the zonal 
velocity component along 185° E. The most interesting 
feature is the secondary maximum occurring near 30°N 
at about 1000 m. This structure, which is located in 
the vicinity of the Hawaiian Ridge becomes a prediction 
which could be tested against in situ observations. Oth- 
erwise, the corrected OGCM shows a large-scale cur- 
rent response with vertically coherent structures over 
the bulk of the section, with the exception of the tropi- 
cal current fluctuations, which tend to be surface inten- 
sified. 

The seasonal current and pressure anomalies are nec- 
essarily associated with temperature (density) changes 
at all depths, some of which are shown in Figure 25 
along the same meridional section at 185° E. Such re- 
sults are equivalent to an estimate of the change in heat 
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Figure 22. Meridional sections of temperature along (a) 165°E and (b) 215°E showing the 
temperature field before (dashed line) and after correction (solid line) by the T/P data. Contour 
increment is 5°C. 


content of the ocean through the seasons. Temperature 
changes of close to 1°C are found in the surface layer. 
Those values are fairly small and should not be con- 
fused with observed surface temperature fluctuations, 
since they represent vertical integrals over the top 100 
m. A simple calculation leads to an associated steric 
height change in sea level 


= -Po 1 



ct66dz 


(28) 


by = 2 cm, consistent with the estimated am- 


plitudes. On the other hand, seasonal temperature 
changes of the order of 0.1° K at 1000 m depth and 
below are artificially large. The relatively low verti- 
cal OGCM resolution enhances the vertical extent of 
responses to surface forcing and exaggerates tempera- 
ture fluctuations in the main thermocline. (Recall the 
vertical resolution discussion in the twin experiments.) 

We claim, then, that the model/data combination 
produces a better estimate of the seasonal cycle than 
either can do alone. The reader may, however, object 
that models such as the GFDL OGCM used here, are 
notoriously unable to properly compute the annual cy- 
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Figure 23. Estimates of seasonal surface elevation anomalies relative to the 1-year mean and 
related geostrophic currents. Fields represent (a) spring, (b), summer, (c) fall, and (d) winter, 
with spring starting at the beginning of March. Positive and negative values are drawn by bold, 
and thin lines, respectively. Contour increment is 1 cm. The reference vector represent 4 cm/s. 
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Figure 24. Meridional section of estimated zonal component of seasonal current anomalies along 
185°E. Fields represent (a) spring, (b), summer, (c) fall, and (d) winter. Contour increment is 
0.2 cm/s in all panels, and negative values are shaded. 


cle of heating and cooling [ Sarmiento , 1986, Stammer ei 
al , 1996]. They generally fail to import sufficient heat 
into the oceanic interior across the sea surface during 
the heating season, thus producing a steric contribution 
to sea level which is too small. In the present circum- 
stances, the observed sea level change is being forced 
onto the OGCM by the Green’s functions. However, 
to the extent that a pressure-excited Green’s function 
mimics the response from thermal driving once the ther- 
mal anomaly has been injected to depth, the results here 
should nonetheless be quite accurate. 

An issue will arise, however, in a model with much 
higher near-surface resolution. Direct thermal injection 
would be expected to penetrate only to about 300 m, 
whereas a mechanically driven perturbation, such as a 
wind-field anomaly would likely (but not necessarily) 
penetrate much deeper. In such a model, Green’s func- 
tion responses for thermal forcing must be computed 
for perturbations with a much shallower vertical pen- 
etration than from mechanical forcing. In the present 


case, the very thick top layer precludes distinguishing 
the two vertical scales, but we believe the results shown 
are as accurate a depiction of the thermal response as 
can be achieved with such low resolution. 

6.3. Residuals 

The estimation method has proceeded on the assump- 
tion that errors in the linear model and data have a 
white noise character. In general this is not correct, 
and residuals, although consistent with the general er- 
ror budget, show strong geographically correlated struc- 
tures, and which at various locations greatly exceed the 
assumed error level. I n Figure 26 the time-averaged 
residuals n = y 0 bs — Y ran ^ ^ = ^00 are shown. 
Values of up to 10 cm in amplitude can be found in var- 
ious locations, but predominantly in the tropics. These 
residuals are static elements of the data which are not 
consistent with the underlying error assumptions and 
linear model physics. 

Such geographically correlated components in the 
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Figure 25. Same as Figure 24, but for estimated seasonal temperature fluctuation. Contour 
increment is 0.1°C in all panels, and negative values are shaded. 


residuals can be caused by various processes. The spa- 
tial patterns actually seen, however, are strongly sug- 
gestive of systematic errors in the T/P data, i.e. geoid 
errors. Stammer and Wunsch [1994] proposed the pres- 
ence of strong geoid errors in the western tropical Pa- 
cific, and the inability of the model to reduce the residu- 
als seen in Figure 26 in that area is supporting evidence 
for the hypothesis. The major competing hypothesis 
would be that the linear model dynamics are inadequate 
to reproduce the observations. Although the analysis is 
not reproduced here, crude estimates suggest that even 
with this low spatial resolution, the model error is in- 
adequate to explain the residuals. If this conclusion is 
accepted, the residuals depicted in Figure 26 can be 
subtracted from the existing geoid estimate to produce 
an improved geoid. This procedure can be quantified 
and formalized and involves using the regionally vary- 
ing estimates of the correlated geoid errors as well as 
the production of a proper error estimate for the model 
(see the discussion of Wunsch and Gaposchkin , [1980], 
on the geoid/circulation estimation problem). 


7. Summary and Discussion 

We have explored the estimation of the ocean cir- 
culation from combined surface altimeter observations, 
a general circulation model, and an associated linear 
model defined from pressure Green’s functions. The 
Green’s function method was applied to two different 
types of data in the Pacific Ocean. Synthetic data from 
an identical twin experiment confirmed the potential of 
the method to estimate the subsurface state of the ocean 
from the sea surface measurements alone. Reasonable 
results were also found for state variables which are not 
directly observed by the altimetry: velocity and temper- 
ature fields over the full water column. For observation 
durations of 1-year or less, there is a loss in resolution 
in the deepest parts of the water column, a result an- 
ticipated by Munk and Wunsch [1982], who proposed 
that acoustic tomography would be the natural com- 
plementary technology. The inference of the need for 
tomography remains a valid one. What is less clear is 
whether the deep ocean will become “observable” by 
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Figure 26. Time mean residuals h = 6p — 6 p Q in 
equivalent surface elevation obtained at K — 500 and 
plotted in geographical order. Contour increment is 2 
cm. 

altimetry alone if the data span can be extended to 
decades, when deep ocean structures would have time 
to manifest themselves at the sea surface. 

The main focus here has been on estimating the 
state of the Pacific Ocean as obtained from 1-year of 
TOPEX/POSEIDON data, with emphasis on the time 
average and the annual cycle. Corrections are found 
to the OGCM initial estimates for these components 
of the circulation, which are both statistically accept- 
able relative to prior estimates, and which also pass the 
test of physical reasonableness. Residuals lying outside 
the prior statistical uncertainties are strongly geograph- 
ically correlated and are almost surely a reflection of 
geoid errors. 

At the surface, the resulting time average state re- 
flects a smoothed version of the T/P data, with reduced 
initial dynamically inconsistent anomalies in the data. 
The correction has weakened the subtropical gyre in 
the model and confined it mainly to the western por- 
tion of the basin. The resulting velocity and tempera- 
ture fields show horizontal and vertical structures which 
are consistent with a generally weakened (cooled) sub- 
tropical gyre, and the change in mean sea surface is 
necessarily accompanied by a change in the mean strat- 
ification. Despite the shortcomings of the OGCM in 
computing the annual heating cycle, the seasonal cycle 
present in the T/P observations is successfully imposed 
on the model, leading to sea surface height changes sim- 
ilar in amplitude and pattern to those observed by T/P, 
with a low, cold North Pacific during winter and spring, 
and a warm, high sea surface during summer and fall. 
Associated current fluctuations are small, except in the 
tropical area, where most of the seasonal changes are in- 
duced not by local heating, but by changing wind fields. 
Seasonal fluctuations in surface elevation of the order of 
a few centimeters are also associated with temperature 
fluctuations of about 0.5°, which show realistic patterns 
and amplitudes. 

These results are preliminary and can be improved in 
a number of ways. Use of information on the geograph- 
ically correlated geoid errors should lead to a significant 
overall reduction in error. Model spatial resolution and 
that of the corresponding Green’s functions should be 


greatly increased. Moreover, we have not taken full ad- 
vantage of all the information available. For example, 
one should use the twice-daily wind and buoyancy flux 
fields available from meteorological analyses to improve 
the OGCM first estimate. Finally, this first attempt 
was entirely focussed on altimetric data, but there is no 
reason not to use the same method with all other data 
types available in the ocean at the same time. In partic- 
ular, the use of deep constraint information as available 
from floats and acoustic tomography should produce 
major improvements in the estimates at depth. 
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